Setup

library(SeqArray)
library(SNPRelate)
library(pander)
library(scales)
library(magrittr)
library(tidyverse)
library(readxl)
library(sp)
library(ggmap)
library(rgdal)
library(ggsn)
library(parallel)
library(qqman)
library(ggrepel)
library(plyranges)
library(rtracklayer)
library(AnnotationDbi)
library(GO.db)
theme_set(theme_bw())
panderOptions("missing", "")
mc <- min(12, detectCores() - 1)
alpha <- 0.05
maxKb <- 40
w <- 169/25.4
h <- 169/25.4

This workflow immediately follows 03_SNPFiltering and performs an analysis on the SNPs retained after filtering. Analysis performed is:

gdsPath <- file.path("..", "5_stacks", "gds", "populations.snps.gds")
gdsFile <- seqOpen(gdsPath, readonly = TRUE)
keepSNPs <- readRDS("keepSNPsAfterLDPruning.RDS")
sampleID <- tibble(
    Sample = seqGetData(gdsFile, "sample.annotation/Sample"),
    Population = seqGetData(gdsFile, "sample.annotation/Population"),
    Location = seqGetData(gdsFile, "sample.annotation/Location")
)
popSizes <- sampleID %>% 
    group_by(Population) %>%
    summarise(n = dplyr::n())
snps <- tibble(
    variant.id = seqGetData(gdsFile, "variant.id") ,
    chromosome = seqGetData(gdsFile, "chromosome"),
    position = seqGetData(gdsFile, "position"),
    snpID = seqGetData(gdsFile, "annotation/id")
) %>%
    mutate(
      snpID = str_remove(snpID, ":.$"),
      snpID = str_replace(snpID, ":", "_")
    ) %>%
    separate(snpID, into = c("Locus ID", "Col"), remove = FALSE) %>%
    mutate_at(vars(`Locus ID`, Col), as.integer) %>%
    mutate(
      Col = Col - 1,
      snpID = paste(`Locus ID`, Col, sep = "_")
    ) %>%
  right_join(keepSNPs)
seqSetFilter(gdsFile, variant.id = snps$variant.id)
## # of selected variants: 18,886
genotypes <- snps %>%
    cbind(
      seqGetData(gdsFile, "genotype") %>% 
        apply(MARGIN = 3, colSums) %>%
        t %>%
        set_colnames(sampleID$Sample)) %>%
  as_tibble() %>%
  gather(Sample, Genotype, one_of(sampleID$Sample)) %>%
  dplyr::filter(!is.na(Genotype)) %>%
  arrange(variant.id, Sample) %>%
  left_join(sampleID)
seqResetFilter(gdsFile)
## # of selected samples: 146
## # of selected variants: 146,763
snpIn1996 <- genotypes %>% 
    filter(Population == 1996) %>% 
    group_by(variant.id) %>% 
    summarise(maf = mean(Genotype) / 2) %>%
    filter(maf > 0)
genotypes %<>% 
    right_join(snpIn1996)
seqClose(gdsFile)

In addition to the above, the set of genes corresponding the Ensembl Release 96 were loaded.

ensGenes <- file.path("..", "external", "Oryctolagus_cuniculus.OryCun2.0.96.gff3.gz") %>%
    import.gff3(feature.type = "gene", sequenceRegionsAsSeqinfo = TRUE) %>%
    .[,c("gene_id", "Name", "description")]

The set of SNPs under investigation was also defined as a GRanges object.

snpsGR <- makeGRangesFromDataFrame(
    df = snps, 
    keep.extra.columns = TRUE, 
    ignore.strand = TRUE, 
    seqinfo = seqinfo(ensGenes), 
    seqnames.field = "chromosome", 
    start.field = "position", 
    end.field = "position"
)

This was intersected with the set of genes to connect SNPs to genes within 40 kb.

snp2Gene <- snpsGR %>% 
    resize(width = 1000*maxKb/2 - 1) %>%
    trim() %>%
    join_overlap_inner(ensGenes)

A set of SNPs within exonic regions was also defined.

ensExons <- file.path("..", "external", "Oryctolagus_cuniculus.OryCun2.0.96.gff3.gz") %>%
    import.gff3(feature.type = "exon", sequenceRegionsAsSeqinfo = TRUE) %>%
  join_overlap_inner_within_directed(
    ensGenes, maxgap = 0, minoverlap = 0, suffix = c(".exon", "")
  ) %>%
  join_overlap_inner(snpsGR)
ensExons <- ensExons[,setdiff(colnames(mcols(ensExons)), "Name.exon")]

PCA

lowCall <- c("gc2901", "gc2776", "gc2731", "gc2727", "gc2686")
snp4PCA <- genotypes %>% 
    filter(!Sample %in% lowCall) %>%
    group_by(variant.id, Population) %>% 
    summarise(n = dplyr::n()) %>% 
    spread(Population, n) %>%
    mutate(N = `1996` + `2010` + `2012`) %>% 
    ungroup() %>%
    filter(N > 0.95*(sum(popSizes$n) - length(lowCall)))
pca <- genotypes %>%
    filter(variant.id %in% snp4PCA$variant.id) %>%
    dplyr::select(variant.id, Sample, Genotype) %>%
    spread(Sample, Genotype) %>%
    as.data.frame() %>%
    column_to_rownames("variant.id") %>%
    as.matrix() %>%
    apply(2, function(x){
        x[is.na(x)] <- mean(x, na.rm = TRUE)
        x
    }) %>%
    t() %>%
    .[, apply(., 2, function(x){length(unique(x)) > 1})] %>%
    prcomp( center = TRUE)

As noted in the previous section sample samples gc2901, gc2776, gc2731, gc2727 and gc2686 had a SNP identification rate \(< 50\)% and as such these were marked as potential outliers. Ignoring these samples, and restricting data to SNPs identified in \(>95\)% of all samples, a preliminary PCA was performed This amounted to 7,767 of the possible 18,878 SNPs for analysis using PCA. Missing values were specified as the mean MAF over all populations combined.

Given the initially observed structure, in which samples from the 2012 population are separating from the other samples which group with the 1996 population, the collection points for the 2012 samples as investigated.

sampleID %<>%
    left_join(
        file.path("..", "external", "GPS_Locations.xlsx") %>%
            read_excel() %>%
            dplyr::select(Sample, ends_with("tude")) %>%
            mutate(Sample = gsub("[Oo][Rr][Aa] ([0-9ABC]*)", "ora\\1", Sample))
    )
*PCA showing population structure. Point size reflects the proportion of SNPs for which imputation was required, and the observed structure appeared independent of this.*

PCA showing population structure. Point size reflects the proportion of SNPs for which imputation was required, and the observed structure appeared independent of this. 

loc <- c(
  range(sampleID$Longitude, na.rm = TRUE) %>% mean,
  range(sampleID$Latitude, na.rm = TRUE) %>% mean
)
saPoly <- readRDS(file.path("..", "external", "saPoly.RDS"))
roads <- readRDS(file.path("..", "external", "roads.RDS"))
gc <- SpatialPoints(cbind(x = 138.655972, y = -31.200305))
proj4string(gc) <- "+proj=longlat +ellps=GRS80 +no_defs"
xBreaks <- seq(138.65, 138.8, by = 0.05)
xLabs <- parse(text = paste(xBreaks, "*degree ~ E", sep = ""))
yBreaks <- seq(-31.2, -31.32, by = -0.04)
yLabs <- parse(text = paste(-yBreaks, "*degree ~ S", sep = ""))
leftN <- tibble(
  x = c(138.7965, 138.8, 138.8) - 0.01,
  y = c(-31.2, -31.198, -31.193)
)
rightN <- tibble(
  x = c(138.8, 138.8, 138.8035) - 0.01,
  y = c(-31.193, -31.198, -31.2)
)
ggMap <- get_map(loc, zoom = 12, maptype = "terrain", color = "bw")
zoomPlot <- ggmap(ggMap, extent = "normal", maprange = FALSE) +
    geom_path(
        aes(long, lat, group = group), 
        data = subset(roads, SURFACE == "UNSE"), 
        linetype = 2, size = 0.3) + 
    geom_path(
        aes(long, lat, group = group), 
        data = subset(roads, SURFACE != "UNSE"), 
        linetype = 1, size = 0.4) +
    geom_label(x = 138.74, y = -31.29, label = "Flinders Ranges NP", alpha = 0.4) +
    geom_label(x = 138.72, y = -31.22, label = "Gum Creek", alpha = 0.4) +
    geom_point(aes(x, y), data = as.data.frame(gc), shape = 3, size = 3) +
    geom_text(aes(x, y), label = "HS", data = as.data.frame(gc), nudge_y = 0.005) +
    geom_polygon(
        data = subset(saPoly, F_CODE == "HD"),
        aes(long, lat, group = group),
        fill = rgb(1, 1, 1, 0), colour = "grey10", size = 0.3) +
    geom_polygon(
        data = subset(saPoly, F_CODE == "PARK"),
        aes(long, lat, group = group),
        fill = rgb(1, 1, 1, 0), colour = "grey10", size = 0.2) +
  geom_point(
    data= filter(pca4Plot, grepl("2012", Population)),
    aes(Longitude, Latitude, colour = Population),
    size = 0.9*ps) +
  geom_polygon(
    aes(x, y), data = leftN, fill = "white", colour = "grey10", size = 0.4
  ) +
  geom_polygon(
    aes(x, y), data = rightN, fill = "grey10", colour = "grey10", size = 0.4
  ) +
  geom_text(
    x = 138.79, y = -31.19, label = "N", colour = "grey10", size = 4
  ) +
  scale_colour_manual(values = popCols[2:3]) +
  coord_cartesian(
    xlim = c(138.618, 138.81),
    ylim = c(-31.335, -31.18),
    expand = 0
  ) +
  scale_x_continuous(breaks = xBreaks, labels = xLabs) +
  scale_y_continuous(breaks = yBreaks, labels = yLabs) +
  guides(colour = FALSE) +
  ggsn::scalebar(
    x.min = 138.618,
    x.max = 138.81,
    y.min = -31.335, 
    y.max = -31.18,
    transform = TRUE,
    dist = 2, 
    dist_unit = "km",
    model = 'GRS80',
    height = 0.012, 
    st.size = 4,
    location = 'bottomright',
    anchor = c(x = 138.8, y = -31.328)) +
  labs(x = "Longitude", y = "Latitude") +
  theme(
    text = element_text(size = fs),
    plot.margin = unit(c(1, 1, 1, 1), "mm")
  )
ausPolygon <- readRDS(file.path("..", "external", "ausPolygon.RDS"))
ausPts <- SpatialPoints(cbind(x = loc[1], y = loc[2]))
proj4string(ausPts) <- proj4string(ausPolygon)
ausPlot <- ggplot() +
  geom_polygon(
    data = ausPolygon, 
    aes(long, lat, group = group),
    fill = "white", colour = "black"
  ) + 
  geom_point(data = as.data.frame(ausPts), aes(x, y), size = 1.5) +
  theme_void() +
  theme(plot.background = element_rect(fill = "white", colour = "black"))
## png 
##   2
*Figure 1: Collection points for all 2012 samples with colours showing sub-populations initially defined by PCA analysis and *k*-means clustering.*

Figure 1: Collection points for all 2012 samples with colours showing sub-populations initially defined by PCA analysis and k-means clustering. 

zoomLoc <- c(138.753, -31.242)
xBreaks <- seq(138.74, 138.76, by = 0.01)
xLabs <- parse(text = paste(xBreaks, "*degree ~ W", sep = ""))
yBreaks <- seq(-31.235, -31.25, by = -0.005)
yLabs <- parse(text = paste(-yBreaks, "*degree ~ S", sep = ""))
central <- rbind(
  x = c(138.749, 138.755),
  y = c(-31.2365, -31.2495)
) %>%
  set_colnames(c("min", "max"))
leftN <- tibble(
  x = c(138.7595, 138.76, 138.76),
  y = c(-31.234, -31.2335, -31.2325)
)
rightN <- tibble(
  x = c(138.76, 138.76, 138.7605),
  y = c(-31.2325, -31.2335, -31.234)
)
get_map(zoomLoc, zoom = 15, maptype = "terrain", source = "google", color = "bw") %>%
  ggmap() +
  geom_jitter(
    data = filter(pca4Plot, grepl("2012", Population)), 
    aes(Longitude, Latitude, colour = Population), 
    size = 3, width = 0.0005, height = 0) +
  geom_rect(
    xmin = central["x", "min"],
    xmax = central["x", "max"],
    ymin = central["y", "min"],
    ymax = central["y", "max"],
    fill = "red", alpha = 0.01, colour = "black") +
  geom_polygon(
    aes(x, y), data = leftN, fill = "white", colour = "grey10", size = 0.4
  ) +
  geom_polygon(
    aes(x, y), data = rightN, fill = "grey10", colour = "grey10", size = 0.4
  ) +
  geom_text(
    x = 138.76, y = -31.232, label = "N", colour = "grey10", size = 5
  ) +
  scale_colour_manual(values = popCols[2:3]) +
  scale_x_continuous(breaks = xBreaks, labels = xLabs) +
  scale_y_continuous(breaks = yBreaks, labels = yLabs) +
  theme_bw() +
  guides(colour = FALSE) +
  labs(x = "Longitude", y = "Latitude") +
  coord_cartesian(
    xlim = c(138.74, 138.762),
    ylim = c(-31.253, -31.231),
    expand = 0
  ) +
  ggsn::scalebar(
    x.min = 138.74, x.max = 138.762, 
    y.min = -31.253, y.max = -31.231,
    transform = TRUE,
    dist = 0.25, dist_unit = "km",, model = 'GRS80', 
    height = 0.012, st.size = 4,
    location = 'bottomright',
    anchor = c(x = 138.761, y = -31.252)
    )
*Zoomed-in view of the central region for 2012 samples with colours showing sub-populations defined by PCA analysis. The region considered to be the Central Region is shaded in red. Due to overlapping GPS points a small amount of jitter has been added to the x-axis.*

Zoomed-in view of the central region for 2012 samples with colours showing sub-populations defined by PCA analysis. The region considered to be the Central Region is shaded in red. Due to overlapping GPS points a small amount of jitter has been added to the x-axis.

Region Analysis

Removal of SNPs Associated with Collection Region

The structure observed within the 2012 population in the PCA could possibly be explained by recent migration into this region. As the samples collected in the outer regions appeared very similar to the 1996 population in the above plots, this would possibly indicate migration a very recent event as the genetic influence of this has not spread through the wider area. Although this may be due to other factors such as sampling bias, this structure was addressed by identifying SNPs which showed an association with the sub-populations identified by PCA analysis. In this way, any candidate SNPs obtained below will be less impacted by this structure, and will be more reflective of the intended variable under study, i.e. selection over time, as opposed to any internal structure of the 2012 population.

oraRegions <- pca4Plot %>% 
  filter(grepl("2012", Population)) %>% 
  rowwise() %>%
  mutate(
    yCentral = cut(Latitude, breaks = central["y",], include.lowest = TRUE), 
    xCentral = cut(Longitude, breaks = central["x",], include.lowest = TRUE),
    Central = (is.na(yCentral) + is.na(xCentral)) == 0
  ) %>%
  dplyr::select(Sample, Central)

Testing for Structure in 2012

This model tests:
H0: No association between genotypes and collection region
HA: An association exists between genotypes and collection region

regionResults <- genotypes %>%
    filter(Population == 2012) %>%
    split(f = .$variant.id) %>%
    mclapply(
      function(x){
        ft <- list(p.value = NA)
        if (length(unique(x$Genotype)) > 1) {
          ft <- x %>%
            left_join(oraRegions) %>%
            group_by(Genotype, Central) %>%
            tally() %>%
            spread(Genotype, n, fill = 0) %>%
            column_to_rownames("Central") %>%
            fisher.test()
        }
        x %>%
          distinct(variant.id, snpID, chromosome, position) %>%
          mutate(p = ft$p.value)
      }, 
      mc.cores = mc
    ) %>%
  bind_rows() %>%
  filter(!is.na(p)) %>%
  arrange(p)

A total of 1682 SNPs were detected as showing a significant association between genotype and the collection region. Under H0, the number expected using α = 0.05 would be 943, and as this number was approximately double that expected, this was taken as evidence of this being a genuine point of concern for this dataset.

Notably, Type II errors were of principle concern in this instance, and as such every SNP with p < 0.05 in the above test was excluded from downstream analysis.

regionSNPs <- filter(regionResults, p < 0.05)
saveRDS(regionSNPs, "regionSNPs.RDS")

Under this additional filtering step, the original set of 18784 SNPs will be reduced to 17,102 for testing by genotype and allele frequency.

Verification Of Removal

In order to verify that the removal of the above SNPs removed the undesired population structure from the 2012 population, the above PCA was repeated, excluding the SNPs marked for removal. The previous structure noted in the data was no longer evident, and as such, these SNPs were marked for removal during analysis by genotype and allele frequency.

pcaPost <- genotypes %>%
  filter(
    variant.id %in% snp4PCA$variant.id,
    !variant.id %in% regionSNPs$variant.id
  ) %>%
  dplyr::select(variant.id, Sample, Genotype) %>%
  spread(Sample, Genotype) %>%
  as.data.frame() %>%
  column_to_rownames("variant.id") %>%
  as.matrix() %>%
  apply(2, function(x){
    x[is.na(x)] <- mean(x, na.rm = TRUE)
    x
  }) %>%
  t() %>%
  .[, apply(., 2, function(x){length(unique(x)) > 1})] %>%
  prcomp( center = TRUE) 
pcaPost4Plot <- pcaPost$x %>%
  as.data.frame() %>%
  rownames_to_column("Sample") %>%
  as_tibble() %>%
  dplyr::select(Sample, PC1, PC2, PC3) %>%
  left_join(sampleID) %>%
  mutate(Cluster = kmeans(cbind(PC1, PC2, PC3), 3)$cluster) %>%
  group_by(Cluster) %>%
  mutate(maxY = max(Latitude, na.rm = TRUE)) %>%
  ungroup() %>%
  mutate(
    Population = case_when(
      Population == 1996 ~ "1996",
      Population == 2010 ~ "Outgroup (Turretfield)",
      maxY == max(maxY) ~ "2012 (Outer)",
      maxY != max(maxY) ~ "2012 (Central)"
    )
  ) %>%
  left_join(
    genotypes %>% 
      filter(variant.id %in% snp4PCA$variant.id, !Sample %in% lowCall) %>% 
      group_by(Sample) %>% 
      tally() %>% 
      mutate(imputationRate = 1 - n / nrow(snp4PCA))
  )
*Figure 2: Principal Components Analysis showing structures before removal of SNPs denoting collection region in the 2012 population, and after removal of these SNPs*

Figure 2: Principal Components Analysis showing structures before removal of SNPs denoting collection region in the 2012 population, and after removal of these SNPs

SNP Analysis

Genotype Frequency Model

This model tests:
H0: No association between genotypes and populations
HA: An association exists between genotypes and populations

genotypeResults <- genotypes %>%
  filter(
    Population != 2010,
    !variant.id %in% regionSNPs$variant.id
  ) %>%
  group_by(variant.id, snpID, Population, Genotype) %>%
  tally() %>%
  ungroup() %>%
  split(f = .$variant.id) %>%
  mclapply(function(x){
    ft <- list(p.value = NA)
    if (length(unique(x$Genotype)) > 1) {
      ft <- x %>%
        spread(Genotype, n, fill = 0) %>%
        column_to_rownames("Population") %>%
        dplyr::select(-variant.id, -snpID) %>%
        fisher.test()
    }
    x %>% 
      distinct(variant.id) %>%
      mutate(p = ft$p.value)
  },
  mc.cores = mc
  ) %>%
  bind_rows() %>%
  filter(!is.na(p)) %>%
  mutate(FDR = p.adjust(p, "fdr"), adjP = p.adjust(p, "bonferroni")) %>%
  arrange(p) %>%
  left_join(
    genotypes %>%
      distinct(variant.id, snpID, chromosome, position)
  ) %>%
  dplyr::select(variant.id, snpID, chromosome, position, everything())
genoSNPs <- genotypeResults %>%
  dplyr::filter(adjP < alpha) %>%
  .[["snpID"]]
SNPs with adjusted p-values < 0.05 when analysing by genotype. Any genes within 40kb of each SNP are indicated.
snpID chromosome position gene_id Name dist2Gene p FDR adjP
3391201_92 GL018705 1,937,359 ENSOCUG00000016428 RCBTB2 37,689 6.103e-08 0.0009035 0.001047
686773_29 3 90,851,512 ENSOCUG00000009616 CRISPLD1 CDS 1.053e-07 0.0009035 0.001807
916731_91 4 89,602,973 ENSOCUG00000006126 RIC8B Intronic 1.865e-07 0.0009837 0.003201
836950_151 4 35,111,855 ENSOCUG00000012144 TMPRSS12 5,745 2.587e-07 0.0009837 0.00444
321186_117 2 27,995,273 2.866e-07 0.0009837 0.004919
2227006_109 13 129,650,867 ENSOCUG00000022823 ID3 16,426 4.419e-07 0.001051 0.007583
3040789_114 19 30,689,426 ENSOCUG00000010440 MSI2 Intronic 4.838e-07 0.001051 0.008303
1089925_112 7 35,584,375 9.805e-07 0.001583 0.01683
464226_128 2 125,603,458 1.014e-06 0.001583 0.01741
3374849_109 GL018704 1,294,009 1.387e-06 0.001984 0.02381
3695415_108 GL018751 737,318 ENSOCUG00000003672 TRAF3 Intronic 1.837e-06 0.002395 0.03153
2099634_108 13 45,086,624 2.028e-06 0.002395 0.0348
1902686_98 12 60,113,885 ENSOCUG00000001871 12,820 2.094e-06 0.002395 0.03593

Under the full genotype model:

  • 15 genotypes were detected as being significantly associated with the two populations when controlling the FWER at α = 0.05
  • 38 genotypes were detected as being significantly associated with the two populations when controlling the FDR at α = 0.05
  • For the most highly ranked SNP (3391201_92), the minor allele has been completely lost in the 2012 population
Figure 3a: Manhattan plot showing results for all SNPs on chromosomes 1:21 for analysis by genotype. The horizontal line indicates the cutoff for an FDR of 5%, with SNPs with an adjusted p-value < 0,05 shown in green

Figure 3a: Manhattan plot showing results for all SNPs on chromosomes 1:21 for analysis by genotype. The horizontal line indicates the cutoff for an FDR of 5%, with SNPs with an adjusted p-value < 0,05 shown in green

*Variants detected as significant to an FDR of 5%. Those denoted with an asterisk received a Bonferroni-adjusted p-value < 0.05 and these typically involved reduction of the minor allele frequency and a shift towards homozygous reference. The remaining sites showed a combination of the same and increasing minor allele abundance, with some variants becoming exclusively heterozygous in the 2012 population.*

Variants detected as significant to an FDR of 5%. Those denoted with an asterisk received a Bonferroni-adjusted p-value < 0.05 and these typically involved reduction of the minor allele frequency and a shift towards homozygous reference. The remaining sites showed a combination of the same and increasing minor allele abundance, with some variants becoming exclusively heterozygous in the 2012 population. 

write_tsv(genotypeResults, "genotypeResults.tsv")

Allele Frequency Model

This model tests:
H0: No association between allele frequencies and populations
HA: An association exists between allele frequencies and populations

alleleResults <- genotypes %>%
  filter(
    Population != 2010,
    !variant.id %in% regionSNPs$variant.id
  ) %>%
  group_by(snpID, Population) %>%
  summarise(
    P = sum(2 - Genotype),
    Q = sum(Genotype)
  ) %>%
  ungroup() %>%
  split(f = .$snpID) %>%
  mclapply(function(x){
    m <- as.matrix(x[c("P", "Q")])
    ft <- list(p.value = NA) 
    if (length(m) == 4) {
      ft <- fisher.test(m)
    }
    x %>%
      mutate(MAF = Q / (P + Q)) %>%
      dplyr::select(snpID, Population, MAF) %>%
      spread(Population, MAF) %>%
      mutate(p = ft$p.value)
  },mc.cores = mc) %>%
  bind_rows() %>%
  filter(!is.na(p)) %>%
  dplyr::rename(
    MAF_1996 = `1996`,
    MAF_2012 = `2012`
  ) %>%
  mutate(
    FDR = p.adjust(p, "fdr"),
    adjP = p.adjust(p, "bonferroni"),
    # Make the MAF is relative to the 1996 population not the reference
    MAF_2012 = case_when(
      MAF_1996 > 0.5 ~ 1 - MAF_2012,
      MAF_1996 <= 0.5 ~ MAF_2012
    ),
    MAF_1996 = case_when(
      MAF_1996 > 0.5 ~ 1 - MAF_1996,
      MAF_1996 <= 0.5 ~ MAF_1996
    )
  ) %>%
  arrange(p) %>%
  left_join(
    genotypes %>% distinct(snpID, chromosome, position)
  ) %>%
  dplyr::select(snpID, chromosome, position, everything())
alleleSnps <- alleleResults %>%
  dplyr::filter(FDR < alpha) %>%
  .[["snpID"]]
SNPs considered as significant when analysing by genotype using an FDR cutoff of 5.00% Minor allele frequencies are reported based on frequencies in the 1996 population. Any genes within 40kb are also shown.
snpID chromosome position MAF_1996 MAF_2012 gene_id Name dist2Gene p FDR adjP
3391201_92 GL018705 1,937,359 0.2255 0.009615 ENSOCUG00000016428 RCBTB2 37,689 3.44e-07 0.005911 0.005916
916731_91 4 89,602,973 0.1961 0 ENSOCUG00000006126 RIC8B Intronic 6.875e-07 0.005911 0.01182
3040789_114 19 30,689,426 0.1863 0 ENSOCUG00000010440 MSI2 Intronic 1.531e-06 0.007628 0.02632
321186_117 2 27,995,273 0.25 0.0122 1.866e-06 0.007628 0.03209
1009106_1 6 9,037,073 0.4432 0.1277 ENSOCUG00000027083 CLEC19A 15,180 2.52e-06 0.007628 0.04333
836950_151 4 35,111,855 0.2907 0.04082 ENSOCUG00000012144 TMPRSS12 5,745 3.086e-06 0.007628 0.05307
464226_128 2 125,603,458 0.1915 0 3.105e-06 0.007628 0.0534
2227006_109 13 129,650,867 0.2708 0.03333 ENSOCUG00000022823 ID3 16,426 4.172e-06 0.008366 0.07174
1902686_98 12 60,113,885 0.1633 0 ENSOCUG00000001871 12,820 4.749e-06 0.008366 0.08166
3374849_109 GL018704 1,294,009 0.2065 0.0102 4.865e-06 0.008366 0.08366
4044777_119 GL018933 204,058 0.01754 0.2019 ENSOCUG00000005859 FNBP1 Intronic 7.749e-06 0.0111 0.1333
1089925_112 7 35,584,375 0.2841 0.04348 9.725e-06 0.01286 0.1672
3695415_108 GL018751 737,318 0.2609 0.03261 ENSOCUG00000003672 TRAF3 Intronic 1.275e-05 0.01447 0.2193
2099634_108 13 45,086,624 0.2609 0.03409 1.412e-05 0.01447 0.2429
87917_119 1 64,635,286 0.186 0 ENSOCUG00000007758 CEP78 Intronic 1.43e-05 0.01447 0.2459
4010189_98 GL018907 283,766 0.1932 0.01 ENSOCUG00000008870 BRWD1 Intronic 1.627e-05 0.01495 0.2798
2689075_49 16 56,527,308 0.05814 0.3043 ENSOCUG00000000984 ESRRG Intronic 1.652e-05 0.01495 0.284
4149102_107 GL019077 34,265 0.2232 0.02941 ENSOCUG00000024268 IFT140 Intronic 2.279e-05 0.01959 0.3919
1997286_93 12 141,844,290 0.08889 0.3556 ENSOCUG00000004829 ESR1 Intronic 2.435e-05 0.01987 0.4187
2906428_79 18 25,311,176 0.1731 0.4434 ENSOCUG00000011533 RHOBTB1 24 2.542e-05 0.01987 0.4371
2405359_66 14 97,813,263 0.02083 0.2021 ENSOCUG00000005835 STXBP5L Intronic 4.48e-05 0.0335 0.7704
3000704_65 19 13,055,154 0.03191 0.234 5.416e-05 0.03776 0.9314
4146664_90 GL019084 74,038 0.1961 0.02041 ENSOCUG00000010071 GNB1 Intronic 5.49e-05 0.03776 0.9441
4098854_101 GL018985 129,495 0.1702 0.01042 ENSOCUG00000006396 Intronic 6.686e-05 0.04259 1
4098838_14 GL018985 123,328 0.3636 0.1275 ENSOCUG00000006396 Intronic 6.797e-05 0.04259 1
1882774_17 12 40,509,779 0.3214 0.5943 7.288e-05 0.04259 1
1045185_75 7 2,702,919 0.2143 0.03774 ENSOCUG00000027861 3,276 7.807e-05 0.04259 1
3016258_99 19 19,178,148 0.22 0.03261 ENSOCUG00000015254 SEZ6 Intronic 7.861e-05 0.04259 1
2028473_73 13 2,979,383 0.01818 0.1731 8.051e-05 0.04259 1
3829050_109 GL018791 957,960 0.163 0.009615 ENSOCUG00000007849 FHAD1 Intronic 8.094e-05 0.04259 1
2580149_42 15 87,473,844 0.3511 0.1154 ENSOCUG00000001910 ADGRL3 Intronic 8.416e-05 0.04259 1
2527760_104 15 43,156,117 0.1531 0.4057 8.594e-05 0.04259 1
3999128_75 GL018883 283,912 0.3654 0.125 ENSOCUG00000016318 SERPINA6 Intronic 8.683e-05 0.04259 1
4176850_23 GL019154 883 0.234 0.04082 ENSOCUG00000002484 MTHFSD 828 8.916e-05 0.04259 1
686773_29 3 90,851,512 0.4375 0.1667 ENSOCUG00000009616 CRISPLD1 CDS 9.731e-05 0.04522 1
3032107_92 19 26,020,372 0.1939 0.02083 ENSOCUG00000024664 LHX1 20,607 0.0001019 0.04612 1

Under this model:

  • 5 SNP alleles were detected as being significantly associated with the two populations when controlling the FWER at α = 0.05. However, as these SNPs were within 21nt of each other, this may represent the same haplotype
  • 38 SNP alleles were detected as being significantly associated with the two populations when controlling the FDR at α = 0.05
## png 
##   2
Manhattan plot showing results for all SNPs on chromosomes 1:21 when analysing by allele frequencies. The horizontal line indicates the cutoff for an FDR of 5%, with SNPs considered significant under the Bonferroni adjustment shown in green.

Manhattan plot showing results for all SNPs on chromosomes 1:21 when analysing by allele frequencies. The horizontal line indicates the cutoff for an FDR of 5%, with SNPs considered significant under the Bonferroni adjustment shown in green.

*Comparison of minor allele frequencies in both populations. SNPs with an FDR-adjusted p-value < 0.05 are coloured red, whilst those with a Bonferroni-adjusted p-value are additionally labelled.*

Comparison of minor allele frequencies in both populations. SNPs with an FDR-adjusted p-value < 0.05 are coloured red, whilst those with a Bonferroni-adjusted p-value are additionally labelled.

*Variants detected as significant to an FDR of 5% under the Allele Frequency model. Those denoted with an asterisk received a Bonferroni-adjusted p-value < 0.05 and these typically involved reduction of the minor allele frequency and a shift towards homozygous reference. The remaining sites showed a combination of the same and increasing minor allele abundance.*

Variants detected as significant to an FDR of 5% under the Allele Frequency model. Those denoted with an asterisk received a Bonferroni-adjusted p-value < 0.05 and these typically involved reduction of the minor allele frequency and a shift towards homozygous reference. The remaining sites showed a combination of the same and increasing minor allele abundance. 

write_tsv(alleleResults, "alleleResults.tsv")

Analysis Using the FLK model

This was performed separately and no SNPs of specific interest were detected.

Export of Data for Bayescan

A VCF was required with only the 1996 and 2012 populations, and restricted to the candidate SNPs after pruning for linkage disequilibrium and detection of the allele in the 1996 population.

var2VCF <- genotypes %>% 
  distinct(variant.id, snpID) %>%
  filter(snpID %in% filter(regionResults, p > 0.05)$snpID) %>% 
  .[["variant.id"]]
gdsFile <- seqOpen(gdsPath, readonly = TRUE)
seqSetFilter(
  gdsFile, 
  variant.id = var2VCF,
  sample.id = sampleID %>% 
    filter(Population %in% c(1996, 2012)) %>% 
    .[["Sample"]]
)
seqGDS2VCF(gdsFile, "../5_stacks/vcf/filtered.vcf.gz")
seqResetFilter(gdsFile)
seqClose(gdsFile)

Inclusion of Bayescan Results

After running Bayescan, results were imported. Unfortunately, these results had IDs which were taken from the line in the above VCF and need to be converted back to IDs compatible with previous analysis.

bayRes <- file.path("..", "external", "BayescanOut.txt.gz") %>%
  gzfile() %>%
  read_delim(delim = " ", col_names = FALSE, skip = 1, col_types = "innnnn-") %>%
  set_colnames(c("ID", "prob","log10.PO.","qval","alpha","fst")) %>%
  mutate(variant.id = var2VCF[ID])
*Summary of results from all three methods, using an FDR of 0.05 for each approach. All SNPs detected under Bayescan were identified using Fisher's Exact Test for both allele and genotype analysis*

Summary of results from all three methods, using an FDR of 0.05 for each approach. All SNPs detected under Bayescan were identified using Fisher’s Exact Test for both allele and genotype analysis 

sigSnps <- c(
  filter(alleleResults, FDR < 0.05)$snpID,
  filter(genotypeResults, FDR < 0.05)$snpID
) %>%
  unique()

Enrichment testing

fullPaths <- list(
  toTable(GOBPANCESTOR) %>% set_names(c("go_id", "ancestor")),
  toTable(GOCCANCESTOR) %>% set_names(c("go_id", "ancestor")),
  toTable(GOMFANCESTOR) %>% set_names(c("go_id", "ancestor"))
) %>% 
  bind_rows() %>%
  dplyr::filter(!ancestor %in% c("all", "GO:0008150", "GO:0003674", "GO:0005575")) %>%
  as_tibble() %>%
  bind_rows(
    tibble(
      go_id = unique(.$go_id),
      ancestor = unique(.$go_id)
    )
  ) %>%
  arrange(go_id, ancestor)
goSummaries <- url("https://uofabioinformaticshub.github.io/summaries2GO/data/goSummaries.RDS") %>%
  readRDS() %>%
  dplyr::rename(go_id = id)

GO Analysis By Gene

dist <- 1000*maxKb

For the initial GO analysis, any gene within 40kb of a SNP considered as significant were tested for any enrichment of biological characteristics.

allSnpsGR <- snps %>%
  dplyr::select(chromosome, position, snpID) %>%
  makeGRangesFromDataFrame(
    keep.extra.columns = TRUE,
    seqnames.field = "chromosome", 
    start.field = "position", 
    end.field = "position", 
    seqinfo = seqinfo(ensGenes)) %>%
  join_overlap_inner(
    ensGenes,
    maxgap = dist
  )
sigSnpGR <- allSnpsGR %>%
  subset(snpID %in% sigSnps)
nGenes <- allSnpsGR$gene_id %>%
  unique() %>%
  length()
nSigGenes <- sigSnpGR$gene_id %>%
  unique() %>%
  length()

10821 genes were identified as being within 40kb of a SNP. For the 32 candidate SNPs within this distance of a gene, 54 genes were identified, and these were analysed for enrichment of any GO terms.

Mappings from genes to GO terms were downloaded manually from the biomart server at www.ensembl.org.

bm <- file.path("..", "external", "gene2GO_biomart.tsv.gz") %>%
  gzfile() %>%
  read_tsv() %>%
  dplyr::filter(!is.na(`GO domain`)) %>%
  dplyr::rename(gene_id = `Gene stable ID`,
         ontology = `GO domain`,
         go_id = `GO term accession`,
         go_term = `GO term name`) %>%
  dplyr::filter(gene_id %in% allSnpsGR$gene_id) %>%
  left_join(fullPaths) %>%
  dplyr::select(gene_id, go_id = ancestor) %>%
  distinct(gene_id, go_id) %>%
  dplyr::filter(!is.na(go_id))

Complete paths of GO terms back to ontology roots were obtained for each gene. The minimum number of steps back to the ontology roots were also obtained for each GO ID.

goRes <- bm %>%
  mutate(inSig = gene_id %in% sigSnpGR$gene_id) %>%
  group_by(go_id) %>%
  summarise(N = dplyr::n(), 
            nSig = sum(inSig)) %>%
  left_join(goSummaries) %>%
  filter(nSig > 0, shortest_path > 2) %>%
  distinct(go_id, .keep_all = TRUE) %>%
  split(f = .$go_id) %>%
  mclapply(function(x){
    ft <- matrix(
      c(nSigGenes - x$nSig, x$nSig, nGenes - nSigGenes - x$N, x$N),
      nrow = 2, byrow = TRUE) %>%
      fisher.test()
    mutate(x, p = ft$p.value)
  }, mc.cores = mc) %>%
  bind_rows() %>%
  mutate(Expected = nSigGenes * N / nGenes,
         FDR = p.adjust(p, "fdr"),
         adjP = p.adjust(p, "bonferroni"),
         term = Term(go_id),
         ontology = Ontology(go_id)) %>%
  dplyr::select(go_id, term, ontology, N, nSig, Expected, p, adjP, FDR) %>%
  arrange(p)
goRes %>%
  dplyr::filter(FDR == min(FDR), nSig > 1) %>%
  mutate(Term = Term(go_id),
         Genes  = vapply(go_id, function(x){
           subset(sigSnpGR, gene_id %in% I(
             dplyr::filter(bm, go_id == x) %>%
               .[["gene_id"]]))$Name %>%
             sort() %>%
             unique() %>%
             paste(collapse = ", ")
         }, character(1))) %>%
  dplyr::select(ID = go_id, Ontology = ontology, Term, Total = N, Sig = nSig, Expected, p, FDR, Genes) %>%
  mutate(ID = gsub(":", "\\\\:", ID),
         Expected = round(Expected, 3),
         p = round(p, 4),
         FDR = round(FDR, 3)) %>%
  distinct(Ontology, Genes, .keep_all = TRUE) %>%
  pander(split.tables = Inf, style = "rmarkdown", 
         justify = "lllrrrrrl",
         caption = paste("Most highly ranked GO terms.", 
                         "Where multiple terms resulted from the identical combination of genes, only the most highly ranked terms is shown.",
                         "Overall, this corresponds to an FDR of", 
                         percent(max(.$FDR))))
Most highly ranked GO terms. Where multiple terms resulted from the identical combination of genes, only the most highly ranked terms is shown. Overall, this corresponds to an FDR of 21.6%
ID Ontology Term Total Sig Expected p FDR Genes
GO:0005496 MF steroid binding 34 3 0.17 8e-04 0.216 ESR1, ESRRG, SERPINA6
GO:0060068 BP vagina development 7 2 0.035 9e-04 0.216 ESR1, LHX1
GO:0001655 BP urogenital system development 146 5 0.729 9e-04 0.216 ESR1, FOXF1, ID3, IFT140, LHX1
GO:0072189 BP ureter development 10 2 0.05 0.0016 0.216 FOXF1, LHX1
GO:0002221 BP pattern recognition receptor signaling pathway 51 3 0.255 0.0024 0.216 ESR1, TRAF3, UFD1
GO:0060041 BP retina development in camera-type eye 56 3 0.279 0.0031 0.216 GNB1, IFT140, LHX1
GO:0043010 BP camera-type eye development 122 4 0.609 0.0035 0.216 COL8A2, GNB1, IFT140, LHX1
GO:0021680 BP cerebellar Purkinje cell layer development 16 2 0.08 0.0036 0.216 LHX1, SEZ6
GO:0072001 BP renal system development 132 4 0.659 0.0046 0.216 FOXF1, ID3, IFT140, LHX1
GO:0003727 MF single-stranded RNA binding 19 2 0.095 0.0048 0.216 AGO3, MSI2
GO:0031050 BP dsRNA processing 19 2 0.095 0.0048 0.216 AGO3, ESR1
GO:0004879 MF nuclear receptor activity 20 2 0.1 0.0053 0.216 ESR1, ESRRG
GO:0050688 BP regulation of defense response to virus 20 2 0.1 0.0053 0.216 TRAF3, UFD1
GO:0001750 CC photoreceptor outer segment 23 2 0.115 0.0068 0.216 GNB1, IFT140
GO:0048701 BP embryonic cranial skeleton morphogenesis 26 2 0.13 0.0085 0.216 IFT140, LHX1
GO:0030522 BP intracellular receptor signaling pathway 85 3 0.424 0.0096 0.216 ESR1, ESRRG, UFD1

Meta Analysis With Previous Paper

The results in the Turretfield population from the paper Resistance to RHD virus in wild Australia rabbits: Comparison of susceptible and resistant individuals using a genomewide approach were then compared to these results. In this previous analysis, results were generated by using Fisher’s Exact Test to test genotypes for association with RHDV susceptibility or resistance.

tfGR <- url("https://raw.githubusercontent.com/hdetering/orycun/master/results/filteredSNPScores.csv") %>% 
  read_csv() %>%
  distinct(min_snp, .keep_all = TRUE) %>%
  dplyr::select(chromosome = chr_name, pos, snpID = min_snp, p, FDR) %>%
  mutate(chromosome = gsub("^chr", "", chromosome)) %>%
  makeGRangesFromDataFrame(
    keep.extra.columns = TRUE,
    ignore.strand = TRUE, 
    seqinfo = seqinfo(ensGenes),
    seqnames.field = "chromosome", 
    start.field = "pos", 
    end.field = "pos") %>%
  sort()

Of the 9,229 unique SNPs contained in the previous dataset, only 1,431 were also identified in the current dataset. A meta-analysis was performed on these SNPs using the results from both separate analyses which used Fisher’s Exact Test on genotypes. In this approach, Fisher’s method was used to combine p-values, and resultant p-values were drawn from a \(chi^2\) distribution with 4 degrees of freedom (i.e. \(2k\)).

combinedResults <- tfGR %>%
  find_overlaps(allSnpsGR) %>%
  mcols() %>%
  as.data.frame() %>%
  as_tibble() %>%
  dplyr::select(snpID = snpID.y, tf.p = p) %>%
  distinct(snpID, .keep_all = TRUE) %>%
  left_join(
    genotypeResults %>%
      dplyr::select(snpID, chromosome, position, geno.p = p)) %>% 
  mutate(
    X = -2*(log(tf.p) + log(geno.p)), 
    p = pchisq(X, 4, lower.tail = FALSE), 
    adjP = p.adjust(p, "bonferroni"),
    FDR = p.adjust(p, "fdr")) %>% 
  arrange(p)

Using Bonferroni’s adjustment to control the FWER at 0.05 gave 2 candidate SNPs, whilst using the more relaxed criterion of an FDR < 0.05 gave 10 candidate SNPs.

Significant SNPs from the meta-analysis using the previous published results from a separate population.
snpID Chromosome Position p adjP FDR
644062_1 3 56309981 1.668e-06 0.002176 0.002176
3378006_1 GL018704 2271135 2.355e-05 0.0307 0.01535
3941015_20 GL018847 216450 5.316e-05 0.06932 0.02311
213788_4 1 146274492 8.646e-05 0.1127 0.02477
3118382_68 20 14614986 9.496e-05 0.1238 0.02477
3108287_38 20 10514947 0.0001406 0.1833 0.03055
3156991_0 21 6185579 0.0001792 0.2337 0.03339
2425370_49 14 114833115 0.0002082 0.2715 0.03394
3904662_77 GL018828 225354 0.0003003 0.3916 0.04325
3970912_82 GL018864 1306 0.0003317 0.4325 0.04325
Genes within 40kb of SNPs identified in the meta-analysis.
snpID Chromosome BP ID Name Description
213788_4 1 146274492 ENSOCUG00000000568 HBB2 hemoglobin, beta
213788_4 1 146274492 ENSOCUG00000021180 HBG2 hemoglobin, gamma G
213788_4 1 146274492 ENSOCUG00000027533 HBG1 hemoglobin, gamma A
213788_4 1 146274492 ENSOCUG00000026912 OR51I1 olfactory receptor family 51 subfamily I member 1
213788_4 1 146274492 ENSOCUG00000025404 OLFR64_1 MOR 5’beta3
644062_1 3 56309981 ENSOCUG00000004827 SFXN1 sideroflexin 1
2425370_49 14 114833115 ENSOCUG00000011289 ALCAM activated leukocyte cell adhesion molecule
3108287_38 20 10514947 ENSOCUG00000008859 KCNH5 potassium voltage-gated channel subfamily H member 5
3118382_68 20 14614986 ENSOCUG00000008913 RPS6KA5 ribosomal protein S6 kinase A5
3156991_0 21 6185579 ENSOCUG00000005961 KDM2B lysine demethylase 2B
3156991_0 21 6185579 ENSOCUG00000005957 RNF34 ring finger protein 34
3378006_1 GL018704 2271135 ENSOCUG00000011042 COL8A2 collagen type VIII alpha 2 chain
3378006_1 GL018704 2271135 ENSOCUG00000011038 ADPRHL2 ADP-ribosylhydrolase like 2
3378006_1 GL018704 2271135 ENSOCUG00000029722 TEKT2 tektin 2
3378006_1 GL018704 2271135 ENSOCUG00000006255 AGO3 argonaute RISC catalytic component 3
3904662_77 GL018828 225354 ENSOCUG00000027620 PAQR4 progestin and adipoQ receptor family member 4
3904662_77 GL018828 225354 ENSOCUG00000017181 PKMYT1 protein kinase, membrane associated tyrosine/threonine 1
3904662_77 GL018828 225354 ENSOCUG00000027985 TNFRSF12A TNF receptor superfamily member 12A
3904662_77 GL018828 225354 ENSOCUG00000017197 CLDN6 claudin 6
3941015_20 GL018847 216450 ENSOCUG00000012298 ODF2 outer dense fiber of sperm tails 2
3941015_20 GL018847 216450 ENSOCUG00000008855 CERCAM cerebral endothelial cell adhesion molecule
3970912_82 GL018864 1306 ENSOCUG00000021601 SCARB1 scavenger receptor class B member 1
Comparison of p-values for SNPs detected in both analyses, with thos considered as significant after meta-analysis shown in red.

Comparison of p-values for SNPs detected in both analyses, with thos considered as significant after meta-analysis shown in red.

Allele distributions in the 1996 and 2012 populations for SNPs identified under the meta-analysis.

Allele distributions in the 1996 and 2012 populations for SNPs identified under the meta-analysis.

GO Enrichment

This list was also used to test for GO enrichment.

sigSnpGR <- allSnpsGR %>%
  subset(snpID %in% filter(combinedResults, FDR < 0.05)$snpID)
nGenes <- allSnpsGR %>%
  subset(snpID %in% combinedResults$snpID) %>%
  unique() %>%
  length()
nSigGenes <- sigSnpGR$gene_id %>%
  unique() %>%
  length()
combinedGoRes <- bm %>%
  dplyr::filter(
    gene_id %in% subset(allSnpsGR, snpID %in% combinedResults$snpID)$gene_id) %>%
  mutate(inSig = gene_id %in% sigSnpGR$gene_id) %>%
  group_by(go_id) %>%
  summarise(N = dplyr::n(), 
            nSig = sum(inSig)) %>%
  left_join(goSummaries) %>%
  filter(nSig > 0, shortest_path > 2) %>%
  distinct(go_id, .keep_all = TRUE) %>%
  split(f = .$go_id) %>%
  mclapply(function(x){
    ft <- matrix(
      c(nSigGenes - x$nSig, x$nSig, nGenes - nSigGenes - x$N, x$N),
      nrow = 2, byrow = TRUE) %>%
      fisher.test()
    mutate(x, p = ft$p.value)
  }, mc.cores = mc) %>%
  bind_rows() %>%
  mutate(Expected = nSigGenes * N / nGenes,
         FDR = p.adjust(p, "fdr"),
         adjP = p.adjust(p, "bonferroni"),
         term = Term(go_id),
         ontology = Ontology(go_id)) %>%
  dplyr::select(go_id, term, ontology, shortest_path, N, nSig, Expected, p, adjP, FDR) %>%
  arrange(p)

Inspection of the results revealed an enrichment for terms related to haemoglobin, however, as these genes are co-located this was revealed to be due to a single SNP with multiple haemoglobin genes within 40kb and can be disregarded.

Most highly ranked GO terms for combined results. Where multiple terms resulted from the identical combination of genes, only the most highly ranked terms is shown. Overall, this corresponds to an FDR of 3.70%
ID Ontology Term Total Sig nSnps Expected p FDR Genes
GO:0015669 BP gas transport 4 3 1 0.075 2e-04 0.034 HBB2, HBG1, HBG2
GO:0019825 MF oxygen binding 4 3 1 0.075 2e-04 0.034 HBB2, HBG1, HBG2
GO:0015893 BP drug transport 14 4 2 0.264 3e-04 0.037 HBB2, HBG1, HBG2, SFXN1

Session Information

R version 3.6.1 (2019-07-05)

Platform: x86_64-pc-linux-gnu (64-bit)

locale: LC_CTYPE=en_AU.UTF-8, LC_NUMERIC=C, LC_TIME=en_AU.UTF-8, LC_COLLATE=en_AU.UTF-8, LC_MONETARY=en_AU.UTF-8, LC_MESSAGES=en_AU.UTF-8, LC_PAPER=en_AU.UTF-8, LC_NAME=C, LC_ADDRESS=C, LC_TELEPHONE=C, LC_MEASUREMENT=en_AU.UTF-8 and LC_IDENTIFICATION=C

attached base packages: stats4, parallel, grid, stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: GO.db(v.3.8.2), AnnotationDbi(v.1.46.1), Biobase(v.2.44.0), rtracklayer(v.1.44.4), plyranges(v.1.4.3), GenomicRanges(v.1.36.1), GenomeInfoDb(v.1.20.0), IRanges(v.2.18.2), S4Vectors(v.0.22.1), BiocGenerics(v.0.30.0), ggrepel(v.0.8.1), qqman(v.0.1.4), ggsn(v.0.5.0), rgdal(v.1.4-4), ggmap(v.3.0.0), sp(v.1.3-1), readxl(v.1.3.1), forcats(v.0.4.0), stringr(v.1.4.0), dplyr(v.0.8.3), purrr(v.0.3.2), readr(v.1.3.1), tidyr(v.0.8.3), tibble(v.2.1.3), ggplot2(v.3.2.1), tidyverse(v.1.2.1), magrittr(v.1.5), scales(v.1.0.0), pander(v.0.6.3), SNPRelate(v.1.18.1), SeqArray(v.1.24.2) and gdsfmt(v.1.20.0)

loaded via a namespace (and not attached): colorspace(v.1.4-1), rjson(v.0.2.20), class(v.7.3-15), futile.logger(v.1.4.3), XVector(v.0.24.0), rstudioapi(v.0.10), bit64(v.0.9-7), lubridate(v.1.7.4), xml2(v.1.2.2), knitr(v.1.24), zeallot(v.0.1.0), jsonlite(v.1.6), Rsamtools(v.2.0.0), broom(v.0.5.2), png(v.0.1-7), compiler(v.3.6.1), httr(v.1.4.1), backports(v.1.1.4), assertthat(v.0.2.1), Matrix(v.1.2-17), lazyeval(v.0.2.2), cli(v.1.1.0), formatR(v.1.7), htmltools(v.0.3.6), tools(v.3.6.1), gtable(v.0.3.0), glue(v.1.3.1), GenomeInfoDbData(v.1.2.1), Rcpp(v.1.0.2), cellranger(v.1.1.0), vctrs(v.0.2.0), Biostrings(v.2.52.0), nlme(v.3.1-141), xfun(v.0.9), rvest(v.0.3.4), XML(v.3.98-1.20), zlibbioc(v.1.30.0), hms(v.0.5.1), SummarizedExperiment(v.1.14.1), lambda.r(v.1.2.3), yaml(v.2.2.0), memoise(v.1.1.0), calibrate(v.1.7.2), stringi(v.1.4.3), RSQLite(v.2.1.2), highr(v.0.8), maptools(v.0.9-5), e1071(v.1.7-2), BiocParallel(v.1.18.1), RgoogleMaps(v.1.4.4), rlang(v.0.4.0), pkgconfig(v.2.0.2), bitops(v.1.0-6), matrixStats(v.0.55.0), evaluate(v.0.14), lattice(v.0.20-38), sf(v.0.7-7), labeling(v.0.3), GenomicAlignments(v.1.20.1), bit(v.1.1-14), tidyselect(v.0.2.5), plyr(v.1.8.4), R6(v.2.4.0), generics(v.0.0.2), DelayedArray(v.0.10.0), DBI(v.1.0.0), pillar(v.1.4.2), haven(v.2.1.1), foreign(v.0.8-72), withr(v.2.1.2), units(v.0.6-4), RCurl(v.1.95-4.12), modelr(v.0.1.5), crayon(v.1.3.4), futile.options(v.1.0.1), KernSmooth(v.2.23-15), rmarkdown(v.1.15), jpeg(v.0.1-8), blob(v.1.2.0), digest(v.0.6.20), classInt(v.0.4-1), VennDiagram(v.1.6.20) and munsell(v.0.5.0)

## [1] TRUE